This exercise was modified based on https://geodacenter.github.io/aot-workshop/Part2-AOT.html.
Goals for this lab: * Deepen the understanding of surface interpolation and the procedure of ordinary kriging. * Practice using R’s variogram and kriging functions (more at: https://gsp.humboldt.edu/olm/R/04_01_Variograms.html) to perform ordinary kriging for temperature point data.
Go through the steps below and understand the meaning of each script line. Then, answer the question at the end of this file.
library(lubridate) #data wrangling
library(sf)
library(sp) #spatial data wrangling & analysis
library(raster) #spatial raster data wrangling
library(terra) #spatial raster data wrangling
#install.packages("gstat")
library(gstat) #kriging and geostatistics
library(tmap) #modern data visualizations
library(leaflet) #modern data visualizations
# Define your working directory
setwd("/Users/jonathanbernard/Desktop/UIUC/Spring 2024/GGIS 224/GGIS224/Labs/Lab11")
nodes <- read.csv("Data_Lab11/AOT_node_temps.txt") #Read the file as csv
head(nodes) #Inspect the data: 31 points with temperature values in both Celsius and Fahrenheit degrees.
## node_id avg_temp avg_temp_f project_id vsn
## 1 001e06109416 26.01265 78.82278 AoT_Chicago 090B
## 2 001e0610b9e9 25.87948 78.58306 AoT_Chicago 80
## 3 001e0610ba8f 28.30287 82.94517 AoT_Chicago 00D
## 4 001e0610bbf9 25.96714 78.74086 AoT_Chicago 20
## 5 001e0610bbff 27.93934 82.29081 AoT_Chicago 25
## 6 001e0610bc10 25.51712 77.93082 AoT_Chicago 01F
## address lat lon description
## 1 Elston and Cortland Chicago IL 41.91659 -87.66641 AoT Chicago (S) [C]
## 2 Broadway Ave & Lawrence Ave Chicago IL 41.96904 -87.65967 AoT Chicago (S) [C]
## 3 Cornell & 47th St Chicago IL 41.81034 -87.59023 AoT Chicago (S)
## 4 Western Ave & 69th St Chicago IL 41.76832 -87.68340 AoT Chicago (S) [C]
## 5 Western Ave & 18th St Chicago IL 41.85780 -87.68581 AoT Chicago (S)
## 6 State St & 87th Chicago IL 41.73631 -87.62418 AoT Chicago (S) [C]
## start_timestamp
## 1 1/1/18 0:00
## 2 2/14/18 0:00
## 3 8/8/17 0:00
## 4 2/13/18 0:00
## 5 12/15/17 0:00
## 6 2/22/18 0:00
Convert the completed node data to spatial object format for plotting and spatial analytics.
node.temps <- nodes #Duplicate 'nodes' to another data file
coordinates(node.temps) <- node.temps[,c("lon", "lat")]
proj4string(node.temps) <- CRS("+init=epsg:4326")
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
There are 31 sensor points with temperature measurements in both Celsius and Fahrenheit degrees.
length(node.temps)
## [1] 31
head(node.temps) #avg_temp for Celsius; avg_temp_f for Fahrenheit.
## node_id avg_temp avg_temp_f project_id vsn
## 1 001e06109416 26.01265 78.82278 AoT_Chicago 090B
## 2 001e0610b9e9 25.87948 78.58306 AoT_Chicago 80
## 3 001e0610ba8f 28.30287 82.94517 AoT_Chicago 00D
## 4 001e0610bbf9 25.96714 78.74086 AoT_Chicago 20
## 5 001e0610bbff 27.93934 82.29081 AoT_Chicago 25
## 6 001e0610bc10 25.51712 77.93082 AoT_Chicago 01F
## address lat lon description
## 1 Elston and Cortland Chicago IL 41.91659 -87.66641 AoT Chicago (S) [C]
## 2 Broadway Ave & Lawrence Ave Chicago IL 41.96904 -87.65967 AoT Chicago (S) [C]
## 3 Cornell & 47th St Chicago IL 41.81034 -87.59023 AoT Chicago (S)
## 4 Western Ave & 69th St Chicago IL 41.76832 -87.68340 AoT Chicago (S) [C]
## 5 Western Ave & 18th St Chicago IL 41.85780 -87.68581 AoT Chicago (S)
## 6 State St & 87th Chicago IL 41.73631 -87.62418 AoT Chicago (S) [C]
## start_timestamp
## 1 1/1/18 0:00
## 2 2/14/18 0:00
## 3 8/8/17 0:00
## 4 2/13/18 0:00
## 5 12/15/17 0:00
## 6 2/22/18 0:00
Confirm the success of spatial object transformation by simple plotting.
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(node.temps) + tm_dots()
We here add the Chicago Community Area spatial data to provide additional context for our maps. The Chicago Community Area shapefile is downloaded from the Chicago Data Portal website: https://data.cityofchicago.org/Facilities-Geographic-Boundaries/Boundaries-Community-Areas-current-/cauq-8yn6.
chiCA <- read_sf('Data_Lab11/ChicagoCommunityAreas/geo_export_fcf7548c-1d07-4c42-99a3-156e428257cc.shp')
Plot the temperature by sensor points, adding Community Areas as a background.
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chiCA) + tm_borders() +
tm_shape(node.temps) + tm_dots(col="avg_temp_f",size=0.3,title="average temp (F)")
We will use (semi)variograms to model the distribution of temperature data. We will then generate a grid on top of the Chicago area, and with the appropriate variogram model selected, use kriging to predict a temperature surface.
A variogram is a function that describes the degree of spatial autocorrelation across data. The final model uses the measure of variability between points at various distances. Points nearby are likely to have more similar values, and as distance between points increases, there is less likely to be similar values between points. In this application, we assume that temperature measurements that are further apart will vary more than measurements taken close together (so, a positive autocorrelation).
The variogram clearly has an outlier, though it may not influence our final predicted surface because the variogram will be fitted by a model.
tmp.vgm <- variogram(node.temps$avg_temp_f ~ 1, node.temps) #here variable avg_temp_f (temperature measurements in Fahrenheit) is used.
plot(tmp.vgm)
We will generate a theoretical spherical model to approximate (fit) the experimental semivariogram.
tmp.fit.sph<- fit.variogram(tmp.vgm, model=vgm("Sph"))
plot(tmp.vgm, tmp.fit.sph)
Next, we’ll create a grid from the Chicago area. The following function will generate a n-by-n grid from a provided spatial data frame.
pt2grid <- function(ptframe,n) {
bb <- bbox(ptframe)
ptcrs <- proj4string(ptframe)
xrange <- abs(bb[1,1] - bb[1,2])
yrange <- abs(bb[2,1] - bb[2,2])
cs <- c(xrange/n,yrange/n)
cc <- bb[,1] + (cs/2)
dc <- c(n,n)
x1 <- GridTopology(cellcentre.offset=cc,cellsize=cs,cells.dim=dc)
x2 <- SpatialGrid(grid=x1,proj4string=CRS(ptcrs))
return(x2)
}
Let’s first generate a grid of 30 by 30 cells using the Chicago Community Area extent. Plot the grid for exploration.
chiCA = as_Spatial(chiCA) #converting chiCA (sf) to a spatial data frame object (needed for executing the pt2grid function)
chi.grid <- pt2grid(chiCA,30)
plot(chi.grid)
To get an even finer resolution, we generate a finer resolution of grid of 100 by 100 cells. Let’s use the finer grid for the following analysis.
chi.grid <- pt2grid(chiCA,100)
plot(chi.grid)
First, we make sure that all our data is in the same projection.
projection(chi.grid) <- crs("+init=epsg:4326")
projection(node.temps) <- crs("+init=epsg:4326")
projection(chiCA) <- crs("+init=epsg:4326")
Krige the data using the spherical model.
temp.kriged <- krige(node.temps$avg_temp_f ~ 1, node.temps, chi.grid, model = tmp.fit.sph) #Again, variable avg_temp_f is used for kriging.
## [using ordinary kriging]
plot(temp.kriged)
Convert the files to SpatRaster and SpatVector files, clip to Chicago boundaries, and plot to confirm everything looks good.
chiCA.spat <- terra::vect(chiCA) #convert so we can use in Terra
temp.kriged.spat <- terra::rast(temp.kriged) #convert so we can use in Terra
chi.temp.kriged <- terra::crop(temp.kriged.spat, chiCA.spat, mask=TRUE) #use chiCA as a mask so we only see temperatures within the Chicago boundary
plot(chi.temp.kriged$var1.pred)
plot(node.temps, add = TRUE)
plot(chiCA, add = TRUE)
Map the kriged Chicago-area surface.
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chiCA) + tm_borders() +
tm_shape(node.temps) + tm_dots(size=0.01) +
tm_shape(chi.temp.kriged$var1.pred) + tm_raster("var1.pred", style = "jenks", title = "Temperature (F)", palette = "BuPu") +
tm_legend(position = c("left", "bottom")) #var1.pred is the ordinary kriging predictions, and #var1.var is the ordinary kriging variance
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
Now it is your turn. Try to krige the temperature values (this time, in Celsius degree) using a exponential fitting model, and map both of the predicted surface and its variance. Hint: to fit the semivariogram by an exponential model, change the argument of model in fit.variogram to “Exp”. Also see more: https://search.r-project.org/CRAN/refmans/gstat/html/fit.variogram.html.
# semivariogram fitting using the exponential model
tmp.vgm <- variogram(node.temps$avg_temp ~ 1, node.temps)
tmp.fit.exp <- fit.variogram(tmp.vgm, model=vgm("Exp"))
plot(tmp.vgm, tmp.fit.exp)
# generating the grid
chi.grid <- pt2grid(chiCA,100)
# preparing the data for kriging analysis
projection(chi.grid) <- crs("+init=epsg:4326")
projection(node.temps) <- crs("+init=epsg:4326")
projection(chiCA) <- crs("+init=epsg:4326")
# krige the data using the exponential model
temp.kriged <- krige(node.temps$avg_temp ~ 1, node.temps, chi.grid, model = tmp.fit.exp)
## [using ordinary kriging]
# Converting the data to to SpatRaster and SpatVector
chiCA.spat <- terra::vect(chiCA)
temp.kriged.spat <- terra::rast(temp.kriged)
chi.temp.kriged <- terra::crop(temp.kriged.spat, chiCA.spat, mask=TRUE)
# mapping the kriged predicted surface using the view mode
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chiCA) + tm_borders() +
tm_shape(node.temps) + tm_dots(size=0.01) +
tm_shape(chi.temp.kriged$var1.pred) + tm_raster("var1.pred", style = "jenks",
title = "Temperature (C)", palette = "BuPu") +
tm_legend(position = c("left", "bottom"))
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Use "fisher" instead of "jenks" for larger data sets
# mapping the kriging variance also using the view mode
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chiCA) + tm_borders() +
tm_shape(node.temps) + tm_dots(size=0.01) +
tm_shape(chi.temp.kriged$var1.var) + tm_raster("var1.var", style = "jenks",
title = "Kriging Variance", palette = "Reds") +
tm_legend(position = c("left", "bottom"))
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Use "fisher" instead of "jenks" for larger data sets
.